Import necessary packages

library(TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(fpp)
## Loading required package: forecast
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA
## Loading required package: fma
## Loading required package: expsmooth
## Loading required package: lmtest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: tseries
library(forecast)
library(tseries)
library(ggplot2)
library(urca)
library(readxl)
library(corrplot)
## corrplot 0.92 loaded
library(stats)

1. Explanatory Analysis

1.1. Import data

df_raw <- read_excel('TS_Project.xlsx')
df <- df_raw

1.2. Check data structure and nulls

head(df)
## # A tibble: 6 × 9
##   time                export import   cpi   ppi  bond sentiment  USDX uncertai…¹
##   <dttm>               <dbl>  <dbl> <dbl> <dbl> <dbl>     <dbl> <dbl>      <dbl>
## 1 2002-01-01 00:00:00   1.41   3.57  178.  128.  5.04      93    120.      117. 
## 2 2002-02-01 00:00:00   1.34   3.52  178   128.  4.91      90.7  119.       87.6
## 3 2002-03-01 00:00:00   1.76   3.79  178.  130.  5.28      95.7  119.       83.4
## 4 2002-04-01 00:00:00   1.46   4.45  179.  131.  5.21      93    115.       88.6
## 5 2002-05-01 00:00:00   1.63   4.36  180.  131.  5.16      96.9  112.       86.3
## 6 2002-06-01 00:00:00   1.76   4.48  180.  131.  4.93      92.4  106.       93.6
## # … with abbreviated variable name ¹​uncertainty
tail(df)
## # A tibble: 6 × 9
##   time                export import   cpi   ppi  bond sentiment  USDX uncertai…¹
##   <dttm>               <dbl>  <dbl> <dbl> <dbl> <dbl>     <dbl> <dbl>      <dbl>
## 1 2022-09-01 00:00:00   6.72   28.5  297.  268.  3.52      58.6  112.       174.
## 2 2022-10-01 00:00:00   5.91   28.4  298.  265.  3.98      59.9  112.       177.
## 3 2022-11-01 00:00:00   5.50   25.6  299.  263.  3.89      56.8  106.       172.
## 4 2022-12-01 00:00:00   5.85   24.8  299.  258.  3.62      59.7  104.       136.
## 5 2023-01-01 00:00:00   6.16   25.0  301.  260.  3.53      64.9  102.       143.
## 6 2023-02-01 00:00:00   5.86   20.6  302.  259.  3.75      67    105.       124.
## # … with abbreviated variable name ¹​uncertainty
any(is.na(df))
## [1] FALSE

1.3. Transform as Time Series data

export = ts(df$export, start = c(2002, 1), frequency = 12)
import = ts(df$import, start = c(2002, 1), frequency = 12)
cpi = ts(df$cpi, start=c(2002,1),frequency=12)
ppi = ts(df$ppi, start=c(2002,1),frequency=12)
bond = ts(df$bond, start=c(2002,1),frequency=12)
sentiment = ts(df$sentiment, start=c(2002,1),frequency=12)
usdx = ts(df$USDX, start=c(2002,1),frequency=12)
uncertainty = ts(df$uncertainty, start=c(2002,1),frequency=12)

1.4. Correlation

# create a correlation matrix
cor_matrix <- cor(cbind(export, import, 
                        cpi, ppi,
                        bond, sentiment,
                        usdx, uncertainty))

# Plot a correlation heatmap
corrplot(cor_matrix, method = "color",
         order = "hclust",
         tl.col = "black", tl.srt = 45, addCoef.col = "black")

# Print the correlation table
# CPI & PPI highly correlated, only keep one.
print(cor_matrix)
##                 export     import        cpi         ppi        bond
## export       1.0000000  0.9576931  0.9508205  0.88258986 -0.61290296
## import       0.9576931  1.0000000  0.9623065  0.91101757 -0.64355565
## cpi          0.9508205  0.9623065  1.0000000  0.92833707 -0.69498684
## ppi          0.8825899  0.9110176  0.9283371  1.00000000 -0.56578495
## bond        -0.6129030 -0.6435556 -0.6949868 -0.56578495  1.00000000
## sentiment   -0.1779160 -0.2435011 -0.2186952 -0.43137804  0.07714775
## usdx         0.2673713  0.2693136  0.2465156  0.01599541 -0.07202055
## uncertainty  0.3651635  0.4453334  0.4491599  0.38277243 -0.59587323
##               sentiment        usdx uncertainty
## export      -0.17791601  0.26737131  0.36516347
## import      -0.24350112  0.26931356  0.44533338
## cpi         -0.21869519  0.24651557  0.44915994
## ppi         -0.43137804  0.01599541  0.38277243
## bond         0.07714775 -0.07202055 -0.59587323
## sentiment    1.00000000  0.32526534 -0.46620569
## usdx         0.32526534  1.00000000  0.04882897
## uncertainty -0.46620569  0.04882897  1.00000000

1.5. Time Series plots for each variable

autoplot(import, main="Import") + ylab('Import(in billion $)')

autoplot(cpi, main="CPI") + ylab('CPI')

autoplot(ppi, main="PPI") + ylab('PPI')

autoplot(bond, main="10Y Govn't Bond") + ylab('Bond')

autoplot(sentiment, main='Consumer Sentiment') + ylab('Sentiment')

autoplot(usdx, main='USD Index') + ylab('USD Index')

autoplot(uncertainty, main='Economic Uncertainty Index') + ylab('Economic Uncertainty Index')

1.6. Stationarity

1.6.1. BoxCox transformation

lambda <- BoxCox.lambda(import)
lambda
## [1] -0.08103269
cbind('Before Transformation' = import, 
      'After Transformation'= BoxCox(import, lambda = lambda)) %>% autoplot(facet = TRUE) + ggtitle('BoxCox Transformation Comparison') + ylab('Import(B$)')

1.6.2. Differencing & KPSS Test

# KPSS Test with the original data
kpss.test(import)
## Warning in kpss.test(import): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  import
## KPSS Level = 3.7518, Truncation lag parameter = 5, p-value = 0.01
# KPSS Test with 1st order differencing data
diff1 = diff(import)
kpss.test(diff1)
## Warning in kpss.test(diff1): p-value greater than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff1
## KPSS Level = 0.027579, Truncation lag parameter = 5, p-value = 0.1
plot(diff1 , type="l", xlab="Year", ylab="1st Diff of GDP", main="Import: 1st Differencing")

# KPSS Test with 1st order differncing + baxcox data
bc_diff1 = diff(BoxCox(import, lambda = lambda))

kpss.test(bc_diff1)
## Warning in kpss.test(bc_diff1): p-value greater than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  bc_diff1
## KPSS Level = 0.10007, Truncation lag parameter = 5, p-value = 0.1
plot(bc_diff1 , type="l", xlab="Year", ylab="1st Diff of GDP", main="Import: 1st Differencing with BoxCox Transformation")

1.6.3. Differencing & ACF, PACF Plot

# Original Data
Acf(import)

Pacf(import)

# Data with 1st order differencing
Acf(diff1)

Pacf(diff1)

# Data with 1st order differencig + boxcox 
Acf(bc_diff1)

Pacf(bc_diff1)

1.7. Check seasonlaity

  • We perform seasonlaity insepctio after making our data stationary, to remove other components to spot seasonlaity in our data more clearly.
# Original Plot
autoplot(import)

# Decomposition Plot (STL)
ts_data <- ts(bc_diff1, frequency = 12)
stl_fit <- stl(ts_data, s.window = "periodic")
plot(stl_fit)

remainder <- stl_fit$time.series[, "remainder"]
plot(remainder)

# Seasonal Plot
ggseasonplot(bc_diff1, polar = TRUE)

ggseasonplot(bc_diff1, polar = FALSE)

# Subseries Plot
ggsubseriesplot(bc_diff1) +
  labs(y = "Import(B$)",
       title = "Seasonal plot: Import")

# Periodogram
spectrum(bc_diff1)

periodogram(bc_diff1)
temp <- periodogram(bc_diff1)

max_freq <- temp$freq[which.max(temp$spec)]
seasonality <- 1/max_freq

1 / temp$freq[43] # 5.953488
## [1] 5.953488
1 / temp$freq[107] #2.392523
## [1] 2.392523
1 / temp$freq[21] # 12.19048
## [1] 12.19048
Acf(bc_diff1)

Pacf(bc_diff1)

2. Model Exploration: ETS

Before COVID 19

2.1. Data preprocessing

2.1.1. Import data

df_raw <- read_excel('TS_Project.xlsx')
df <- ts(df_raw$import, start=c(2002,01), frequency=12)

2.1.2. Plot data

autoplot(df) + ylab('Import') + ggtitle('US Monthly Import(2002/01 ~ 2023/02)')

2.1.3. Train-Test Split

df_train <- window(df, start=c(2002, 1), end=c(2018, 12))
df_test <- window(df, start=c(2019, 1), end=c(2019, 12))

2.2. ETS Modeling

2.2.0 Seasonal Naive(Baseline Model)

fcst_snaive <- snaive(df_train, h = 12)
autoplot(snaive(df_train, h = 12))

accuracy(fcst_snaive)
##                     ME     RMSE      MAE    MPE     MAPE MASE      ACF1
## Training set 0.7617909 1.527044 1.280903 7.1367 12.59435    1 0.7750122
accuracy(fcst_snaive$mean, df_test)
##                 ME     RMSE     MAE      MPE     MAPE       ACF1 Theil's U
## Test set 0.4499564 1.248044 1.03557 2.415875 5.877225 -0.3960617 0.7297898
checkresiduals(fcst_snaive) 

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 682.64, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

2.2.1 Simple Exponential Smoothing

  • Simple exponential smoothing is a time series forecasting method that uses a weighted average of past observations to make future predictions.
  • The method is based on the assumption that the future values of a time series can be estimated as a weighted average of the past values, with more recent observations being given greater weight.
ses3 <- ets(df_train, model='ANN')
fcst_ses3 <- forecast(ses3, h=12)
accuracy(fcst_ses3)[,2] 
## [1] 0.9509481
accuracy(fcst_ses3$mean, df_test)[,2] 
## [1] 1.643461
checkresiduals(fcst_ses3$residuals)

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 172.5, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
autoplot(ses3) + ylab('Import($B)')

autoplot(ses3) + autolayer(fitted(ses3)) + ylab('Import($B)')

2.2.2. Holt’s Exponential Smoothing

  • The ETS model assumes that the observed time series can be decomposed into three components: trend, seasonal, and error.
  • The trend component refers to the long-term pattern or direction of the time series. The seasonal component refers to the repeating pattern that occurs at fixed intervals of time, such as daily, weekly, or yearly. The error component refers to the random and unpredictable variation in the data that cannot be explained by the trend or seasonal components.
# 1) Linear trend with additive seasonality
hwa <- hw(df_train, seasonal = 'additive', h=12, damped=FALSE)
accuracy(hwa) 
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.003863464 0.6488222 0.4852083 -0.2996543 4.800245 0.3788017
##                    ACF1
## Training set 0.01406249
accuracy(hwa$mean, df_test) 
##                 ME      RMSE       MAE       MPE     MAPE       ACF1 Theil's U
## Test set 0.1271767 0.8469332 0.6890742 0.4701921 3.986005 -0.4069091 0.4849101
checkresiduals(hwa$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 49.132, df = 24, p-value = 0.001823
## 
## Model df: 0.   Total lags used: 24
# pdf("plot.pdf", width = 10, height = 6)
plot(hwa)
lines(df_test, col = "red")

# 2) Linear trend with multiplicative seasonality
hwm <- hw(df_train, seasonal = 'multiplicative', h=12, damped=FALSE)
accuracy(hwm) 
##                      ME      RMSE       MAE         MPE     MAPE      MASE
## Training set 0.02226477 0.5818108 0.4363622 -0.01005463 4.344781 0.3406676
##                     ACF1
## Training set -0.03132167
accuracy(hwm$mean, df_test) 
##                 ME      RMSE       MAE       MPE     MAPE       ACF1 Theil's U
## Test set 0.2187908 0.8660047 0.7331691 0.9624011 4.256856 -0.3067592 0.5090764
checkresiduals(hwm$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 42.431, df = 24, p-value = 0.01154
## 
## Model df: 0.   Total lags used: 24
plot(hwm)
lines(df_test, col = "red")

# 3) Linear trend with additive seasonality and damping
hwad <- hw(df_train, seasonal = 'additive', h=12, damped=TRUE)
accuracy(hwad) 
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.08943375 0.6597823 0.4971829 0.8109482 5.319108 0.3881503
##                     ACF1
## Training set -0.01316333
accuracy(hwad$mean, df_test) 
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set 0.2356866 1.105793 0.913462 0.7802561 5.422637 0.1745505 0.6544695
checkresiduals(hwad$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 40.299, df = 24, p-value = 0.01986
## 
## Model df: 0.   Total lags used: 24
plot(hwad)
lines(df_test, col = "red")

# 4) Linear trend with multiplicative seasonality and damping
hwmd <- hw(df_train, seasonal = 'multiplicative', h=12, damped=TRUE)
accuracy(hwmd) 
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.03776369 0.5807925 0.4377261 0.2958402 4.469252 0.3417324
##                     ACF1
## Training set -0.02054866
accuracy(hwmd$mean, df_test) 
##                 ME      RMSE       MAE      MPE     MAPE       ACF1 Theil's U
## Test set 0.4950831 0.9971667 0.8049686 2.566052 4.561863 -0.2181502 0.5892602
checkresiduals(hwmd$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 35.775, df = 24, p-value = 0.0577
## 
## Model df: 0.   Total lags used: 24
plot(hwmd)
lines(df_test, col = "red")

2.2.3. Holt-Winter’s Exponenetial Smoothing

# 1 ETS_BEST
ets_best <- ets(df_train) 

summary(ets_best) 
## ETS(M,A,M) 
## 
## Call:
##  ets(y = df_train) 
## 
##   Smoothing parameters:
##     alpha = 0.5353 
##     beta  = 1e-04 
##     gamma = 2e-04 
## 
##   Initial states:
##     l = 4.0986 
##     b = 0.0424 
##     s = 0.9908 1.0682 1.1106 1.0167 1.0183 1.0399
##            1.0286 0.9756 0.9869 0.9801 0.8526 0.9316
## 
##   sigma:  0.057
## 
##      AIC     AICc      BIC 
## 847.9923 851.2826 904.4003 
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.03831786 0.5906786 0.4378753 0.1381647 4.319084 0.3418489
##                     ACF1
## Training set -0.03521243
checkresiduals(ets_best$residuals)

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 42.295, df = 24, p-value = 0.01196
## 
## Model df: 0.   Total lags used: 24
fcst_ets_best <- forecast(ets_best, h=12)
accuracy(fcst_ets_best$mean, df_test)[,2] #RMSE:2.534726
## [1] 0.8525379
autoplot(forecast(ets_best)) + ylab('Import($B)')

plot(ets_best)

plot(fcst_ets_best)
lines(df_test, col = "red")

# 2. ETS_AAA

ets1 <- ets(df_train, model='AAA')

summary(ets1) 
## ETS(A,A,A) 
## 
## Call:
##  ets(y = df_train, model = "AAA") 
## 
##   Smoothing parameters:
##     alpha = 0.4797 
##     beta  = 1e-04 
##     gamma = 0.2997 
## 
##   Initial states:
##     l = 4.0752 
##     b = 0.0591 
##     s = -0.2185 0.1891 0.3965 -0.0095 0.2708 0.3481
##            0.3534 0.1724 0.1472 -0.1551 -0.8346 -0.6597
## 
##   sigma:  0.6759
## 
##      AIC     AICc      BIC 
## 942.3971 945.6874 998.8051 
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.003863464 0.6488222 0.4852083 -0.2996543 4.800245 0.3788017
##                    ACF1
## Training set 0.01406249
checkresiduals(ets1$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 49.132, df = 24, p-value = 0.001823
## 
## Model df: 0.   Total lags used: 24
fcst_ets1 <- forecast(ets1, h=12)
accuracy(fcst_ets1$mean, df_test)[,2] 
## [1] 0.8469332
autoplot(forecast(ets1)) + ylab('Import($B)')

plot(ets1)

plot(fcst_ets1)
lines(df_test, col = "red")

# 3. ETS MAA
ets2 <- ets(df_train, model='MAA') 

summary(ets2) 
## ETS(M,A,A) 
## 
## Call:
##  ets(y = df_train, model = "MAA") 
## 
##   Smoothing parameters:
##     alpha = 0.5215 
##     beta  = 1e-04 
##     gamma = 0.206 
## 
##   Initial states:
##     l = 4.1148 
##     b = 0.0575 
##     s = 0.0136 0.3469 0.5217 0.1403 0.2663 0.2594
##            0.1876 -0.0763 -0.0638 -0.1533 -0.8209 -0.6217
## 
##   sigma:  0.0615
## 
##      AIC     AICc      BIC 
## 880.9684 884.2588 937.3765 
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.007855664 0.6566577 0.4792403 -0.2734595 4.698261 0.3741425
##                     ACF1
## Training set 0.004567237
checkresiduals(ets2$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 38.875, df = 24, p-value = 0.02815
## 
## Model df: 0.   Total lags used: 24
fcst_ets2 <- forecast(ets2, h=12)
accuracy(fcst_ets2$mean, df_test)[,2] 
## [1] 0.8176148
autoplot(forecast(ets2)) + ylab('Import($B)')

plot(ets2)

plot(fcst_ets2)
lines(df_test, col = "red")

# 4. ETS MMM
ets3 <- ets(df_train, model='MMM') 
summary(ets3) 
## ETS(M,M,M) 
## 
## Call:
##  ets(y = df_train, model = "MMM") 
## 
##   Smoothing parameters:
##     alpha = 0.554 
##     beta  = 1e-04 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 4.037 
##     b = 1.0072 
##     s = 0.9947 1.0738 1.1076 1.0177 1.0197 1.0366
##            1.0237 0.9757 0.9852 0.9811 0.8512 0.933
## 
##   sigma:  0.0567
## 
##      AIC     AICc      BIC 
## 847.6910 850.9813 904.0990 
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.01799567 0.5900927 0.4405025 -0.2720982 4.360535 0.3438999
##                     ACF1
## Training set -0.05031675
checkresiduals(ets3) 

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,M,M)
## Q* = 40.548, df = 24, p-value = 0.01866
## 
## Model df: 0.   Total lags used: 24
fcst_ets3 <- forecast(ets3, h=12)
accuracy(fcst_ets3$mean, df_test)[,2] 
## [1] 0.8680253
autoplot(forecast(ets3)) + ylab('Import($B)')

plot(ets3)

plot(fcst_ets3)
lines(df_test, col = "red")

After COVID19

3.1. Data preprocessing

3.1.1. Train-Test Split

df_train2 <- window(df, start=c(2002, 1), end=c(2021, 12))
df_test2 <- window(df, start=c(2022, 1), end=c(2022,12))

3.1.2 Seasonal Naive(Baseline Model)

fcst_snaive <- snaive(df_train2, h = 12)
autoplot(snaive(df_train2, h = 12))

accuracy(fcst_snaive $mean, df_test2)
##                ME     RMSE      MAE      MPE     MAPE      ACF1 Theil's U
## Test set 3.242702 4.445163 4.111537 12.07819 15.53626 0.4679377  1.611265

3.2. ETS Modeling

3.2.1 Simple Exponential Smoothing

ses32 <- ets(df_train2, model='ANN')
fcst_ses32 <- forecast(ses32, h=12)
accuracy(fcst_ses32)[,2] 
## [1] 1.118179
accuracy(fcst_ses32$mean, df_test2)[,2] 
## [1] 2.103524
checkresiduals(fcst_ses32$residuals)

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 211.62, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
autoplot(ses32) + ylab('Import($B)')

autoplot(ses32) + autolayer(fitted(ses32)) + ylab('Import($B)')

3.2.2. Holt’s Exponential Smoothing

# 1) Linear trend with additive seasonality
hwa2 <- hw(df_train2, seasonal = 'additive', h=12, damped=FALSE)
accuracy(hwa2) 
##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.03557122 0.7700349 0.5754799 -0.2529396 5.100559 0.3898795
##                    ACF1
## Training set 0.01213894
accuracy(hwa2$mean, df_test2) 
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -2.142104 3.432513 2.447793 -8.391731 9.549002 0.4576512  1.273803
checkresiduals(hwa2$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 40.881, df = 24, p-value = 0.01716
## 
## Model df: 0.   Total lags used: 24
plot(hwa2)
lines(df_test2, col = "red")

# 2) Linear trend with multiplicative seasonality
hwm2 <- hw(df_train2, seasonal = 'multiplicative', h=12, damped=FALSE)
accuracy(hwm2)
##                      ME      RMSE       MAE         MPE     MAPE      MASE
## Training set 0.04448036 0.7014787 0.5107636 0.005129007 4.569962 0.3460351
##                    ACF1
## Training set 0.06717551
accuracy(hwm2$mean, df_test2) 
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -1.570837 3.174066 2.495035 -6.166118 9.709162 0.4946886  1.172481
checkresiduals(hwm2$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 60.478, df = 24, p-value = 5.467e-05
## 
## Model df: 0.   Total lags used: 24
plot(hwm2)
lines(df_test2, col = "red")

# 3) Linear trend with additive seasonality and damping
hwad2 <- hw(df_train2, seasonal = 'additive', h=12, damped=TRUE)
accuracy(hwad2) 
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.07959986 0.7972801 0.6008245 0.5927588 5.847758 0.4070501
##                    ACF1
## Training set 0.05177696
accuracy(hwad2$mean, df_test2) 
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -3.092347 3.829269 3.108496 -12.07427 12.13198 0.2978517  1.418878
checkresiduals(hwad2$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 58.148, df = 24, p-value = 0.0001161
## 
## Model df: 0.   Total lags used: 24
plot(hwad2)
lines(df_test2, col = "red")

# 4) Linear trend with multiplicative seasonality and damping
hwmd2 <- hw(df_train2, seasonal = 'multiplicative', h=12, damped=TRUE)
accuracy(hwmd2) 
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.05315586 0.6975741 0.5102869 0.1244391 4.548272 0.3457122
##                    ACF1
## Training set 0.02301767
accuracy(hwmd2$mean, df_test2) 
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -2.176211 3.633945 2.760502 -8.450697 10.74315 0.5221199  1.339611
checkresiduals(hwmd2$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 52.418, df = 24, p-value = 0.0006899
## 
## Model df: 0.   Total lags used: 24
plot(hwmd2)
lines(df_test2, col = "red")

3.2.3. Holt-Winter’s Exponenetial Smoothing

# 1. ETS_BEST
ets_best2 <- ets(df_train2) 

summary(ets_best2) 
## ETS(M,A,M) 
## 
## Call:
##  ets(y = df_train2) 
## 
##   Smoothing parameters:
##     alpha = 0.5311 
##     beta  = 1e-04 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 4.0927 
##     b = 0.0597 
##     s = 1.0063 1.0846 1.1101 1.0194 1.0135 1.0378
##            1.0187 0.9762 0.9892 0.9744 0.8362 0.9336
## 
##   sigma:  0.0581
## 
##      AIC     AICc      BIC 
## 1094.530 1097.287 1153.701 
## 
## Training set error measures:
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.06211406 0.7294272 0.5263934 0.0311818 4.498561 0.3566241
##                    ACF1
## Training set 0.01879104
checkresiduals(ets_best2$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 46.18, df = 24, p-value = 0.004212
## 
## Model df: 0.   Total lags used: 24
fcst_ets_best2 <- forecast(ets_best2, h=12)
accuracy(fcst_ets_best2$mean, df_test2)[,2] 
## [1] 2.019115
autoplot(forecast(ets_best2)) + ylab('Import($B)')

plot(ets_best2)

plot(fcst_ets_best2)
lines(df_test2, col = "red")

# 2. ETS_AAA

ets12 <- ets(df_train2, model='AAA')

summary(ets12) 
## ETS(A,A,A) 
## 
## Call:
##  ets(y = df_train2, model = "AAA") 
## 
##   Smoothing parameters:
##     alpha = 0.4867 
##     beta  = 0.0276 
##     gamma = 0.346 
## 
##   Initial states:
##     l = 4.2205 
##     b = 0.1176 
##     s = 0.202 0.0634 0.5544 -0.1505 -0.2504 0.3038
##            0.2454 0.0457 0.3895 0.1847 -0.6942 -0.8938
## 
##   sigma:  0.7971
## 
##      AIC     AICc      BIC 
## 1223.920 1226.677 1283.091 
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.03557122 0.7700349 0.5754799 -0.2529396 5.100559 0.3898795
##                    ACF1
## Training set 0.01213894
checkresiduals(ets12$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 40.881, df = 24, p-value = 0.01716
## 
## Model df: 0.   Total lags used: 24
fcst_ets12 <- forecast(ets12, h=12)
accuracy(fcst_ets12$mean, df_test2)[,2] 
## [1] 3.432513
autoplot(forecast(ets12)) + ylab('Import($B)')

plot(ets12)

plot(fcst_ets12)
lines(df_test2, col = "red")

# 3. ETS MAA
ets22 <- ets(df_train2, model='MAA') 

summary(ets22) 
## ETS(M,A,A) 
## 
## Call:
##  ets(y = df_train2, model = "MAA") 
## 
##   Smoothing parameters:
##     alpha = 0.5565 
##     beta  = 1e-04 
##     gamma = 0.2409 
## 
##   Initial states:
##     l = 4.3895 
##     b = 0.0618 
##     s = -0.042 0.2136 0.6139 0.1723 0.2651 0.2256
##            0.2715 -0.1323 0.0125 -0.2621 -0.7543 -0.5837
## 
##   sigma:  0.062
## 
##      AIC     AICc      BIC 
## 1126.354 1129.110 1185.524 
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE     MAPE      MASE
## Training set 0.05355054 0.773105 0.5608316 -0.1563987 4.828425 0.3799555
##                    ACF1
## Training set 0.01192749
checkresiduals(ets22$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 35.481, df = 24, p-value = 0.06159
## 
## Model df: 0.   Total lags used: 24
fcst_ets22 <- forecast(ets22, h=12)
accuracy(fcst_ets22$mean, df_test2)[,2] 
## [1] 1.907392
autoplot(forecast(ets22)) + ylab('Import($B)')

plot(ets22)

plot(fcst_ets22)
lines(df_test2, col = "red")

# 4. ETS MMM
ets32 <- ets(df_train2, model='MMM') 
summary(ets32) 
## ETS(M,M,M) 
## 
## Call:
##  ets(y = df_train2, model = "MMM") 
## 
##   Smoothing parameters:
##     alpha = 0.5646 
##     beta  = 1e-04 
##     gamma = 9e-04 
## 
##   Initial states:
##     l = 4.0127 
##     b = 1.0074 
##     s = 1.0144 1.0841 1.1101 1.0237 1.0135 1.0326
##            1.0172 0.9665 0.9824 0.9872 0.838 0.9304
## 
##   sigma:  0.058
## 
##      AIC     AICc      BIC 
## 1094.544 1097.301 1153.715 
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE     MAPE      MASE
## Training set 0.01542948 0.7196395 0.5242459 -0.1484539 4.482695 0.3551692
##                     ACF1
## Training set -0.03362053
checkresiduals(ets32$residuals) 

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 46.808, df = 24, p-value = 0.003535
## 
## Model df: 0.   Total lags used: 24
fcst_ets32 <- forecast(ets32, h=12)
accuracy(fcst_ets32$mean, df_test2)[,2] 
## [1] 2.680248
autoplot(forecast(ets32)) + ylab('Import($B)')

plot(ets32)

plot(fcst_ets32)
lines(df_test2, col = "red")